module evolve

  use dimensions

  complex, private      :: ExpKv(jz, jv)
  complex, private      :: ExpKx(jx, jDir)
  complex, private      :: ExpKy(jy, jDir)
  complex, private      :: ExpKz(jz)

  complex, private      :: uxMod(jx)
  complex, private      :: uyMod(jy)

  real, private         ::  weight 

contains

!-------------------------------------------------------------------------------------
subroutine evolve_init()
!-------------------------------------------------------------------------------------

  use parameters, only : dt, kz, vPhase, ux0, uy0, Dxxxx, Dzzzz, &
                         iuModulation, iunit, Lx, Ly, SideLoss

  use fft_wrap, only   : fft_forward, fft_backward
  use mpi_wrap

  implicit none

  integer              :: iDir, ikz, iv


  !-- setting modulation vectors --

  write(msg,*) char(10), ' Initializing evolve module'
  call mpi_wrap_msg


  if (iuModulation) then

     uxMod(:) = exp( - 14 * xvec(:)**6 / (Lx/2)**6 )
     uyMod(:) = exp( - 14 * yvec(:)**6 / (Ly/2)**6 )

     call fft_forward(jx,uxMod,uxMod)
     uxMod=uxMod/float(jx)
     call deAlias(jx,uxMod)        
     call fft_backward(jx,uxMod,uxMod)

     call fft_forward(jy,uyMod,uyMod)
     uyMod=uyMod/float(jy)
     call deAlias(jy,uyMod)        
     call fft_backward(jy,uyMod,uyMod)

  else

     uxMod(:) = 1.
     uyMod(:) = 1.

  endif

  if (Lx .eq. 0) uxMod(:) = 1.
  if (Ly .eq. 0) uyMod(:) = 1.
 



  !-- initialization --

  ! weight = 0.5 since each direction has 1/2 density,
  ! goes to 1/4 is 3D (1/6 with hexagonal symmetry)

  weight = 1./(jDir*npDir)

  ExpKz(:) = cexp( -iunit * kZvec(:) * (-vPhase) * dt )
 
  do iv=1,jv
     ExpKv(:,iv) = cexp( -iunit*kZvec(:)*vvec(iv)*dt - Dzzzz * kZvec(:)**4 *dt  - SideLoss*dt  )
  enddo


  do iDir = 1,jDir
     ExpKx(:,iDir) = cexp( - iunit * kXvec(:)*ux0(iDir) * dt  )
     ExpKy(:,iDir) = cexp( - iunit * kYvec(:)*uy0(iDir) * dt  )
     !ExpKx(:,iDir) = cexp( - iunit * kXvec(:)*ux0(iDir) * dt - Dxxxx * kXvec(:)**4 *dt )
     !ExpKy(:,iDir) = cexp( - iunit * kYvec(:)*uy0(iDir) * dt - Dxxxx * kYvec(:)**4 *dt )
  enddo

end subroutine evolve_init


!-------------------------------------------------------------------------------------
subroutine evolve_F_long_advect(dfu, halfstep)
!-------------------------------------------------------------------------------------

  use fft_wrap, only   : fft_forward, fft_backward

  implicit none

  real,    intent(inout)   :: dfu(jv, jz, jy, jx, jDir)
  logical, intent(in)      :: halfstep
  complex                  :: fz(jz)
  integer                  :: ix, iy, iz, iv, iDir


 !------------------------------------------

  !-- advect in z: df/dt = -v df/dz + Dzzzz d4f/dz4 - SideLoss

  if (halfstep) then

     do iDir = 1,jDir

        do ix=1,jx
           do iy=1, jy
              do iv=1,jv
                 fz(:) = dfu(iv,:,iy,ix,iDir)
                 call fft_forward(jz,fz,fz)
                 fz = fz/float(jz)
                 call dealias(jz,fz)              
                 fz(:) = fz(:)*sqrt(ExpKv(:,iv))
                 call fft_backward(jz,fz,fz)
                 dfu(iv,:,iy,ix,iDir) = fz(:)
              enddo
           end do
        enddo

     enddo

  else

     do iDir = 1,jDir

        do ix=1,jx
           do iy=1, jy
              do iv=1,jv
                 fz(:) = dfu(iv,:,iy,ix,iDir)
                 call fft_forward(jz,fz,fz)
                 fz = fz/float(jz)
                 call dealias(jz,fz)              
                 fz(:) = fz(:)*ExpKv(:,iv)
                 call fft_backward(jz,fz,fz)
                 dfu(iv,:,iy,ix,iDir) = fz(:)
              enddo
           end do
        enddo

     enddo

  endif

end subroutine evolve_F_long_advect




!-------------------------------------------------------------------------------------
subroutine evolve_F_trans_advect(dfu, &                        ! input/output
                           uxMat, uyMat, f0vec)                ! input
!-------------------------------------------------------------------------------------

  use parameters, only : Dxxxx, dt, iBcPerpPeriodic, ux0, uy0, iLinVlasov

  implicit none

  real,    intent(inout)   :: dfu(jv, jz, jy, jx, jDir)
  real,    intent(in)      :: uxMat(jz, jy, jx, jDir)
  real,    intent(in)      :: uyMat(jz, jy, jx, jDir)
  real,    intent(in)      :: f0vec(jv)

  real                     :: axVar(jx+1), fxVec(jx+1)
  real                     :: ayVar(jy+1), fyVec(jy+1)
  integer                  :: ix, iy, iz, iv, iDir
  integer                  :: jjx, jjy, iBcFlag

 !------------------------------------------

  if (iLinVlasov) then
     call evolve_F_trans_advect_lin(dfu, uxMat, uyMat, f0vec)
     return
  endif


  if (iBcPerpPeriodic) then
     jjx = jx
     jjy = jy
     iBcFlag = 1
  else  
     jjx = jx + 1
     jjy = jy + 1
     iBcFlag = 0
  end if


  do iDir = 1,jDir


     !-- advect in x: d(f+f0)/dt = -d/dx [(f+f0) Ux(x,y)] --

     do iy=1,jy
        do iz=1,jz                      
           axVar(1:jx) = uxMat(iz,iy,:,iDir)           
           axVar(jx+1) = ux0(iDir)          
           do iv=1,jv
              fxVec(1:jx) = dfu(iv,iz,iy,:,iDir) + weight * f0vec(iv)
              fxVec(jx+1) = weight * f0vec(iv)
              call UpWind2ndOrder_Variable_a( fxVec,  dt, dx, axVar, jjx, &
                   weight*f0vec(iv),  weight*f0vec(iv), iBcFlag )
              dfu(iv,iz,iy,:,iDir) = fxVec(1:jx)  - weight * f0vec(iv) 
           enddo
        enddo
     end do

     !-- advect in y: d(f+f0)/dt = -d/dy [(f+f0) Uy(x,y)] --

     do ix=1,jx
        do iz=1,jz                      
           ayVar(1:jy) = uyMat(iz,:,ix,iDir)           
           ayVar(jy+1) = uy0(iDir)
           do iv=1,jv
              fyVec(1:jy) = dfu(iv,iz,:,ix,iDir) + weight * f0vec(iv)
              fyVec(jy+1) = weight * f0vec(iv)
              call UpWind2ndOrder_Variable_a( fyVec,  dt, dy, ayVar, jjy, &
                   weight*f0vec(iv),  weight*f0vec(iv), iBcFlag )
              dfu(iv,iz,:,ix,iDir) = fyVec(1:jy)  - weight * f0vec(iv) 
           enddo
        enddo
     end do


   enddo ! iDir

end subroutine evolve_F_trans_advect




!-------------------------------------------------------------------------------------
subroutine evolve_F_trans_advect_lin(dfu, &                    ! input/output
                           uxMat, uyMat, f0vec)                ! input
!-------------------------------------------------------------------------------------

  use fft_wrap, only   : fft_forward, fft_backward

  use parameters, only:  dt, iunit

  implicit none

  real,    intent(inout)   :: dfu(jv, jz, jy, jx, jDir)
  real,    intent(in)      :: uxMat(jz, jy, jx, jDir)
  real,    intent(in)      :: uyMat(jz, jy, jx, jDir)
  real,    intent(in)      :: f0vec(jv)


  complex    :: dummy1Dx(jx)
  complex    :: dummy1Dy(jy)

  integer    :: ix, iy, iz, iv, iDir

   
  !-- advect dfu in y:  d/dt (dfu) = - Uy0 d/dy(dfu) - f0 d/dy(Uy)


  do iDir = 1,jDir
   
     !-- y-advection, first term: d/dt (dfu) = - Uy0 d/dy(dfu)

     do ix=1,jx
        do iz=1,jz

           do iv=1,jv

              dummy1Dy(:) = dfu(iv,iz,:,ix,iDir)
              call fft_forward(jy,dummy1Dy,dummy1Dy)
              dummy1Dy=dummy1Dy/float(jy)
              call dealias(jy,dummy1Dy) 
              dummy1Dy(:) = ExpKy(:,iDir) * dummy1Dy(:)
              call fft_backward(jy,dummy1Dy,dummy1Dy)
              dfu(iv,iz,:,ix,iDir) = dummy1Dy(:)

           end do

        end do
     end do


     !-- x-advection, first term: d/dt (dfu) = - Ux0 d/dx(dfu)

     do iy=1,jy
        do iz=1,jz


           do iv=1,jv

              dummy1Dx(:) = dfu(iv,iz,iy,:,iDir)
              call fft_forward(jx,dummy1Dx,dummy1Dx)
              dummy1Dx=dummy1Dx/float(jx)
              call dealias(jx,dummy1Dx) 
              dummy1Dx(:) = ExpKx(:,iDir) * dummy1Dx(:)
              call fft_backward(jx,dummy1Dx,dummy1Dx)
              dfu(iv,iz,iy,:,iDir) = dummy1Dx(:)

           end do

        end do
     end do



    !-- y-advection, second term:  d/dt (dfu) = - f0 d/dy(Uy)

     do ix=1,jx
        do iz=1,jz

           dummy1Dy(:) = uyMat(iz,:,ix,iDir)
           call fft_forward(jy,dummy1Dy,dummy1Dy)
           dummy1Dy=dummy1Dy/float(jy)
           call dealias(jy,dummy1Dy)

           dummy1Dy(:) = - iunit * kYvec(:) * dummy1Dy(:) * dt

           call fft_backward(jy,dummy1Dy,dummy1Dy)

           do iv=1,jv
              dfu(iv,iz,:,ix,iDir) =  dfu(iv,iz,:,ix,iDir) + dummy1Dy(:)*weight*f0vec(iv)
           end do

        end do
     enddo


     !-- x-advection, second term:  d/dt (dfu) = - f0 d/dx(Ux)


     do iy=1,jy
        do iz=1,jz

           dummy1Dx(:) = uxMat(iz,iy,:,iDir)
           call fft_forward(jx,dummy1Dx,dummy1Dx)
           dummy1Dx=dummy1Dx/float(jx)
           call dealias(jx,dummy1Dx)

           dummy1Dx(:) = - iunit * kXvec(:) * dummy1Dx(:) * dt

           call fft_backward(jx,dummy1Dx,dummy1Dx)

           do iv=1,jv
              dfu(iv,iz,iy,:,iDir) =  dfu(iv,iz,iy,:,iDir) + dummy1Dx(:)*weight*f0vec(iv)
           end do

        end do
     enddo

  end do ! iDir loop


end subroutine evolve_F_trans_advect_lin




!-------------------------------------------------------------------------------------
subroutine evolve_F_accel(dfu, &                                ! input/output
                          Ez2D, f0vec)                          ! input
!-------------------------------------------------------------------------------------

  use fft_wrap, only   : fft_forward, fft_backward
  use parameters, only : dt, Dvvvv, CollisionRate, iunit

  implicit none

  real,    intent(inout)   :: dfu(jv, jz, jy, jx, jDir)
  real,    intent(in)      :: Ez2D(jz,jy,jx)
  real,    intent(in)      :: f0vec(jv)

  complex                  :: f(jv), f0(jv), fa(jv), b(jv), c(jv), e(jv)
  integer                  :: ix, iy, iz, iv, iDir

  real                     :: rc

  rc = CollisionRate

  f0(:)=f0vec(:)*weight      
  call fft_forward(jv,f0,f0)   !  f=FT(dfu)
  f0=f0/float(jv)              !  normalized velocity  FT of equilibriium


  !-- acceleration: df/dt = Ez df/dv

   do iDir = 1,jDir

      do ix=1,jx
         do iy=1,jy

           do iv=1,jv
               fa(iv) = sum(dfu(iv,:,iy,ix,iDir))/jz;
           enddo
   
           call fft_forward(jv,fa,fa)
           fa=fa/float(jv)
 
           do iz=1,jz   !	begin spatial loop
               f(:) = dfu(:,iz,iy,ix,iDir)
               call fft_forward(jv,f,f)
               f = f/float(jv)
               call dealias(jv,f)
               b(:) =  - Dvvvv*kVvec(:)**4  - rc*kVvec(:)**2 - iunit*kVvec(:)*Ez2D(iz,iy,ix)
               e(:) =    cexp(b(:)*dt)
               c(:) = -iunit*kVvec(:)*Ez2D(iz,iy,ix)*f0 + Dvvvv*kVvec(:)**4*fa
               c(:) = c(:)/b(:)
               c(1) = 0
               f(:) = f*e + c * (e - 1)
               call fft_backward(jv,f,f)
               dfu(:,iz,iy,ix,iDir)=f(:)       
            enddo
         enddo
      enddo

   enddo
   

end subroutine evolve_F_accel

!-------------------------------------------------------------------------------------
subroutine evolve_F_erase_nonlin(dfu, &                    ! input/output
                                 fInit)                    ! input
!-------------------------------------------------------------------------------------

  use fft_3d

  use parameters, only:  modeSave

  implicit none

  real,    intent(inout)   :: dfu(jv, jz, jy, jx, jDir)
  complex, intent(in)      :: fInit(jv,jz,jDir)


  complex    :: dummy3D(jz,jy,jx)
  complex    :: fp(jz)
  complex    :: fm(jz)

  integer    :: iv, iDir, px, py, mx, my


  !-- index of the main mode

  px = modeSave + 1
  mx = jx - modeSave + 1
  py = 1
  my = 1


  do iDir = 1,jDir
     do iv=1,jv

        dummy3D(:,:,:) = dfu(iv,:,:,:,iDir)

        call fft_forward_x(dummy3D)
        call fft_forward_y(dummy3D)

        !-- preserve the main mode

        fp(:)=dummy3D(:,py,px)
        fm(:)=dummy3D(:,my,mx)

        dummy3D(:,:,:) = 0;

        !-- restore the zero's mode and x,y-independent equilibrium
 

        dummy3D(:,py,px) = fp(:)
        dummy3D(:,my,mx) = fm(:)
        dummy3D(:,1,1) = fInit(iv,:,iDir)


        call fft_backward_y(dummy3D)
        call fft_backward_x(dummy3D)

        dfu(iv,:,:,:,iDir) = dummy3D(:,:,:)
 
     end do
  end do


 

end subroutine evolve_F_erase_nonlin




!-------------------------------------------------------------------------------------
subroutine evolve_Uz(dfu, f0vec, &                             ! input
                     uzMat)                                    ! output
!-------------------------------------------------------------------------------------

  implicit none

  real,    intent(in)      :: dfu(jv, jz, jy, jx, jDir)
  real,    intent(in)      :: f0vec(jv)
  real,    intent(out)     :: uzMat(jz, jy, jx, jDir)


  real       :: fVsVvec(jv)
  integer    :: ix, iy, iz, iv, iDir


  !-- integrate over "v" to find longtitudinal velocity Uz -- 
       
  do iDir = 1,jDir

     do ix = 1,jx
        do iy = 1,jy
           do iz = 1,jz 
              fVsVvec(:) = dfu(:,iz,iy,ix,iDir) + weight * f0vec(:)   ! weight or simpy 0.5?
              uzMat(iz,iy,ix,iDir) = sum( fVsVvec * vVec ) / sum( fVsVvec )
           enddo
        enddo
     enddo

  enddo
   

end subroutine evolve_Uz


!-------------------------------------------------------------------------------------
subroutine evolve_U_long_advect(uxMat, uyMat, &                       ! input/output
                                uzMat)                                ! input
!-------------------------------------------------------------------------------------

  use fft_wrap, only   : fft_forward, fft_backward

  use parameters, only : dt, iFourier, iWaveFrame


  implicit none

  real,    intent(inout)   :: uxMat(jz, jy, jx, jDir)
  real,    intent(inout)   :: uyMat(jz, jy, jx, jDir)
  real,    intent(in)      :: uzMat(jz, jy, jx, jDir)


  complex    :: dummy1Dz(jz)

  real       :: uxVecZ(jz)
  real       :: uyVecZ(jz)
  real       :: uzVecZ(jz)

  integer    :: ix, iy, iz, iv, iDir

  if( iFourier==1 .and. iWaveFrame==0) return
  
  !-- advect U in z:  dUx/dt = -Uz dUx/dz, dUy/dt = -Uz dUy/dz
    
  do iDir = 1,jDir
   
     !  advect ux in z since <u_z> approximately = -vphase =  V0, added on June 2, 2009
     do ix=1,jx
        do iy=1,jy

           if( iFourier==1 ) then                   ! FT with linearization in uz

              dummy1Dz(:) = uxMat(:,iy,ix,iDir)
              call fft_forward(jz,dummy1Dz,dummy1Dz)
              dummy1Dz=dummy1Dz/float(jz)
              call dealias(jz,dummy1Dz) 
              dummy1Dz(:) = ExpKz(:) * dummy1Dz(:)
              call fft_backward(jz,dummy1Dz,dummy1Dz)
              uxMat(:,iy,ix,iDir) = dummy1Dz(:)

              dummy1Dz(:) = uyMat(:,iy,ix,iDir)
              call fft_forward(jz,dummy1Dz,dummy1Dz)
              dummy1Dz=dummy1Dz/float(jz)
              call dealias(jz,dummy1Dz) 
              dummy1Dz(:) = ExpKz(:) * dummy1Dz(:)
              call fft_backward(jz,dummy1Dz,dummy1Dz)
              uyMat(:,iy,ix,iDir) = dummy1Dz(:)


           else                                     ! upwind or central differences

              uzVecZ(:) = uzMat(:,iy,ix,iDir)
              uxVecZ(:) = uxMat(:,iy,ix,iDir)
              uyVecZ(:) = uyMat(:,iy,ix,iDir)

              call SubUpwind_NonCon_Form( uxVecZ, dt, dz, uzVecZ, jz, &
                   0.0, 0.0, 1, -iFourier, dummy1Dz )

              call SubUpwind_NonCon_Form( uyVecZ, dt, dz, uzVecZ, jz, &
                   0.0, 0.0, 1, -iFourier, dummy1Dz )

              uxMat(:,iy,ix,iDir)= uxVecZ(:)
              uyMat(:,iy,ix,iDir)= uyVecZ(:)
 
           endif

        end do
     enddo

  end do


end subroutine evolve_U_long_advect

!-------------------------------------------------------------------------------------
subroutine evolve_U_trans_advect_lin(uxMat, uyMat)                    ! input/output
!-------------------------------------------------------------------------------------

  use fft_wrap, only   : fft_forward, fft_backward

  use parameters, only : dt


  implicit none

  real,    intent(inout)   :: uxMat(jz, jy, jx, jDir)
  real,    intent(inout)   :: uyMat(jz, jy, jx, jDir)

  complex    :: dummy1Dx(jx)
  complex    :: dummy1Dy(jy)

  integer    :: ix, iy, iz, iv, iDir

   
  !-- advect Ux and Uy in y:  dUx/dt = -Uy0 dUx/dy,   dUy/dt = -Uy0 dUy/dy

  do iDir = 1,jDir
   
     do ix=1,jx
        do iz=1,jz

           dummy1Dy(:) = uxMat(iz,:,ix,iDir)
           call fft_forward(jy,dummy1Dy,dummy1Dy)
           dummy1Dy=dummy1Dy/float(jy)
           call dealias(jy,dummy1Dy) 
           dummy1Dy(:) = ExpKy(:,iDir) * dummy1Dy(:)
           call fft_backward(jy,dummy1Dy,dummy1Dy)
           uxMat(iz,:,ix,iDir) = dummy1Dy(:)

           dummy1Dy(:) = uyMat(iz,:,ix,iDir)
           call fft_forward(jy,dummy1Dy,dummy1Dy)
           dummy1Dy=dummy1Dy/float(jy)
           call dealias(jy,dummy1Dy) 
           dummy1Dy(:) = ExpKy(:,iDir) * dummy1Dy(:)
           call fft_backward(jy,dummy1Dy,dummy1Dy)
           uyMat(iz,:,ix,iDir) = dummy1Dy(:)

        end do
     enddo

  end do


  !-- advect Ux and Uy in x:  dUx/dt = -Ux0 dUx/dy,   dUy/dt = -Ux0 dUy/dy

  do iDir = 1,jDir
   
     do iy=1,jy
        do iz=1,jz

           dummy1Dx(:) = uxMat(iz,iy,:,iDir)
           call fft_forward(jx,dummy1Dx,dummy1Dx)
           dummy1Dx=dummy1Dx/float(jx)
           call dealias(jx,dummy1Dx) 
           dummy1Dx(:) = ExpKx(:,iDir) * dummy1Dx(:)
           call fft_backward(jx,dummy1Dx,dummy1Dx)
           uxMat(iz,iy,:,iDir) = dummy1Dx(:)

           dummy1Dx(:) = uyMat(iz,iy,:,iDir)
           call fft_forward(jx,dummy1Dx,dummy1Dx)
           dummy1Dx=dummy1Dx/float(jx)
           call dealias(jx,dummy1Dx) 
           dummy1Dx(:) = ExpKx(:,iDir) * dummy1Dx(:)
           call fft_backward(jx,dummy1Dx,dummy1Dx)
           uyMat(iz,iy,:,iDir) = dummy1Dx(:)

        end do
     enddo

  end do


end subroutine evolve_U_trans_advect_lin


!-------------------------------------------------------------------------------------
subroutine evolve_U_trans_advect(uxMat, uyMat)                        ! input/output
!-------------------------------------------------------------------------------------

  use parameters, only : ux0, uy0, dt, iBcPerpPeriodic, iLinHydro

  implicit none

  real,    intent(inout)   :: uxMat(jz, jy, jx, jDir)
  real,    intent(inout)   :: uyMat(jz, jy, jx, jDir)

  real       :: dummy1Dy(jy)
  real       :: dummy1Dx(jx)

  real       :: uxVecX(jx+1)
  real       :: uyVecX(jx+1)

  real       :: uxVecY(jy+1)
  real       :: uyVecY(jy+1)

  integer    :: ix, iy, iz, iv, iDir
  integer    :: jjx, jjy, iBcFlag

  if (iLinHydro) then
    call evolve_U_trans_advect_lin(uxMat, uyMat)
    return
  endif

  if (iBcPerpPeriodic) then
     jjx = jx
     jjy = jy
     iBcFlag = 1
  else  
     jjx = jx + 1
     jjy = jy + 1
     iBcFlag = 0
  end if

 

  !-- advect U in y:  dUx/dt = -Uy dUx/dy,   dUy/dt = -Uy dUy/dy


  do iDir = 1,jDir
   
     do ix=1,jx
        do iz=1,jz
             
           uxVecY(1:jy) = uxMat(iz,:,ix,iDir)
           uyVecY(1:jy) = uyMat(iz,:,ix,iDir)
           uxVecY(jy+1) = ux0(iDir)
           uyVecY(jy+1) = uy0(iDir)

           call SubUpwind_NonCon_Form( uxVecY, dt, dy, uyVecY, jjy, &
                ux0(iDir), ux0(iDir), iBcFlag, 1, dummy1Dy )

           call NonLinearAdvection( uyVecY, dt, dy, jjy, &
                uy0(iDir), uy0(iDir), iBcFlag )  

           uxMat(iz,:,ix,iDir) = uxVecY(1:jy)
           uyMat(iz,:,ix,iDir) = uyVecY(1:jy)

        end do
        enddo
        
  end do

 
  !-- advect U in x:  dUx/dt = -Ux dUx/dx,   dUy/dt = -Ux dUy/dx

  do iDir = 1,jDir

     do iy=1,jy 
        do iz=1,jz

           uxVecX(1:jx) = uxMat(iz,iy,:,iDir)
           uyVecX(1:jx) = uyMat(iz,iy,:,iDir)
           uxVecX(jx+1) = ux0(iDir)
           uyVecX(jx+1) = uy0(iDir)


           call SubUpwind_NonCon_Form( uyVecX, dt, dx, uxVecX, jjx, &
                uy0(iDir), uy0(iDir), iBcFlag, 1, dummy1Dx )

           call NonLinearAdvection( uxVecX, dt, dx, jjx, &
                ux0(iDir), ux0(iDir), iBcFlag )  
 
           uxMat(iz,iy,:,iDir) = uxVecX(1:jx)
           uyMat(iz,iy,:,iDir) = uyVecX(1:jx)

        end do
     end do

  end do

end subroutine evolve_U_trans_advect


!-------------------------------------------------------------------------------------
 subroutine evolve_U_accel(uxMat, uyMat, &                     ! input/output
                           Ex2D, Ey2D)                         ! input
!-------------------------------------------------------------------------------------

  use parameters, only : dt

  implicit none

  real,    intent(inout)   :: uxMat(jz, jy, jx, jDir)
  real,    intent(inout)   :: uyMat(jz, jy, jx, jDir)

  real,    intent(in)      :: Ex2D(jz,jy,jx)
  real,    intent(in)      :: Ey2D(jz,jy,jx)

  integer                  :: iDir

  !-- acceleration:  dUx/dt = Ex,  dUy/dt = Ey


  do iDir = 1,jDir
   
     uxMat(:,:,:,iDir) = uxMat(:,:,:,iDir) + Ex2D * dt
     uyMat(:,:,:,iDir) = uyMat(:,:,:,iDir) + Ey2D * dt
         
  enddo
 

end subroutine evolve_U_accel



!-------------------------------------------------------------------------------------
subroutine evolve_U_erase_nonlin(uxMat, uyMat)             ! input/output

  use fft_3d
  use parameters, only:  modeSave, ux0, uy0

  implicit none

  real,    intent(inout)   :: uxMat(jz, jy, jx, jDir)
  real,    intent(inout)   :: uyMat(jz, jy, jx, jDir)

  complex    :: dummy3D(jz,jy,jx)
  complex    :: fp(jz)
  complex    :: fm(jz)

  integer    :: iDir, px, py, mx, my


  !-- index of the main mode

  px = modeSave + 1
  mx = jx - modeSave + 1
  py = 1
  my = 1


  !-- x-velocity

  do iDir = 1,jDir
 
        dummy3D(:,:,:) = uxMat(:,:,:,iDir) - ux0(iDir) 

        call fft_forward_x(dummy3D)
        call fft_forward_y(dummy3D)

        !-- preserve the main mode

        fp(:)=dummy3D(:,py,px)
        fm(:)=dummy3D(:,my,mx)

        dummy3D(:,:,:) = 0;

        !-- restore the main mode
 
        dummy3D(:,py,px) = fp(:)
        dummy3D(:,my,mx) = fm(:)


        call fft_backward_y(dummy3D)
        call fft_backward_x(dummy3D)

        uxMat(:,:,:,iDir) = dummy3D(:,:,:) + ux0(iDir) 

  end do


  !-- y-velocity

  do iDir = 1,jDir
 
        dummy3D(:,:,:) = uyMat(:,:,:,iDir) - uy0(iDir) 

        call fft_forward_x(dummy3D)
        call fft_forward_y(dummy3D)

        !-- preserve the main mode

        fp(:)=dummy3D(:,py,px)
        fm(:)=dummy3D(:,my,mx)

        dummy3D(:,:,:) = 0;

        !-- restore the main mode
 
        dummy3D(:,py,px) = fp(:)
        dummy3D(:,my,mx) = fm(:)


        call fft_backward_y(dummy3D)
        call fft_backward_x(dummy3D)

        uyMat(:,:,:,iDir) = dummy3D(:,:,:) + uy0(iDir) 

  end do


 

end subroutine evolve_U_erase_nonlin


!-------------------------------------------------------------------------------------
 subroutine evolve_U_filter(uxMat, uyMat, &                     ! input/output
                            filter)                             ! input
!-------------------------------------------------------------------------------------

  use fft_3d

  use parameters, only : ux0, uy0

  implicit none


  real,    intent(inout)   :: uxMat(jz, jy, jx, jDir)
  real,    intent(inout)   :: uyMat(jz, jy, jx, jDir)

  real,    intent(in)      :: filter(jz,jy,jx)

  complex    :: du(jz,jy,jx)

  integer    :: ix, iy, iz, iDir


  do iDir = 1,jDir

     !-- filtering x-velocity

     du = uxMat(:,:,:,iDir) - ux0(iDir) 


     do ix=1,jx
        du(:,:,ix) = du(:,:,ix)*uxMod(ix)
     end do

     do iy=1,jy
        du(:,iy,:) = du(:,iy,:)*uyMod(iy)
     end do


     call fft_forward_x(du)
     call fft_forward_y(du)

     du(:,:,:) = du(:,:,:)*filter(:,:,:)

     call fft_backward_x(du)
     call fft_backward_y(du)

     uxMat(:,:,:,iDir) = du +  ux0(iDir)


     !-- filtering y-velocity

     du = uyMat(:,:,:,iDir) - uy0(iDir) 

     do ix=1,jx
        du(:,:,ix) = du(:,:,ix)*uxMod(ix)
     end do

     do iy=1,jy
        du(:,iy,:) = du(:,iy,:)*uyMod(iy)
     end do

     call fft_forward_x(du)
     call fft_forward_y(du)

     du(:,:,:) = du(:,:,:)*filter(:,:,:)

     call fft_backward_x(du)
     call fft_backward_y(du)

     uyMat(:,:,:,iDir) = du +  uy0(iDir)

  enddo


  if (jx .eq. 1) uxMat = 0.
  if (jy .eq. 1) uyMat = 0.

end subroutine evolve_U_filter
 

!-------------------------------------------------------------------------------------

end module evolve
